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Abstract 

First principles based phase diagram calculations were performed for the hexagonal closest packed 
octahedral-interstitial solid solution system aH fOx (aHf[ ]\-xOx\ [ ]=Vacancy; < X < 1/2). 
The cluster expansion method was used to do a ground state analysis, and to calculate the phase 
diagram. The predicted diagram has four ordered ground-states in the range < X < 1/2, but 
one of these, at X=5/12, is predicted to disproportionate at T« 220K. At Xps 1/3 (HfeO) and 

o 

O Xw 1/2 (Hf20), order-parameter vs temperature plots evince a cascade of ordered structures. 
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The group 4 hexagonal closest packed (hep) suboxides MOx (M = Ti, Zr or Hf) all 
exhibit octahedral interstitial ordering of oxygen, O, and vacancies, [ ], in solid solutions of 
the form aM[ ]i-xOx, < X < 1/2). By far, the most studied of these systems is ZrOx, 
because of issues related to the oxidation of Zircalloy cladding on UO2 fuel rods in light-water 
reactors.^^ The hep-based HfOx system has attracted less attentionp^"^ but Hafnium 
alloys are also potential cladding materials; e.g. for long-lived nuclear waste transmutation 
applications in Boiling Water Reactors.^ Also, investigating the chemical systematics of 
all three group 4 suboxides enhances understanding of each binary system. In the ZrOx 
system, long-period superstructure (LPSS) phases were reported^ in samples with Xr^I/3, 
but not predicted in a recent first principles phase diagram (FPPD) calculation^ however, 
in the HfOx FPPD calculations described below, a cascade of related ordered structures 
is predicted at X^l/3 and X«l/2; it is not yet clear if this cascade constitutes a Devil's 
Staircase. 1221231 

I. METHODOLOGY 

A. Total Energy Calculations 

Formation energies, AEf (Fig. [I]) were calculated for fully relaxed hep aHf, HfO 
(hep aHf with all octahedral interstices occupied by O), and 96 aHf[ ]i- n O n supercells 
of intermediate composition. All calculations were performed with the density functional 
theory (DFT) based Vienna ab initio simulation program (VASP, version 4.4.523EH) using 
projector- augmented plane-wave pseudopotentials, and the generalized gradient approxima- 
tion for exchange and correlation energies. Electronic degrees of freedom were optimized 
with a conjugate gradient algorithm, and both cell constant and ionic positions were fully 
relaxed. 

Total energy calculations were converged with respect to k-point meshes by increasing 
the density of k-points for each structure until convergence is achieved. A 500 eV energy 
cutoff was used, in the "high precision" option which guarantees that absolute energies are 
converged to within a few meV/site (a few tenths of a kJ/site of exchangeable species; O, 
I ]). Residual forces were typically 0.02 eV or less. 

Calculated formation energies, AEf, relative to a mechanical mixture of aHf + aHf'O, 



for the 106 aHf[ ]i- n O n supercells are plotted as solid circles in Fig. flj Values of AEf are, 

AEf = (Egtr — E aH f — E aH f )/(2) (1) 

where: E S t r is the total energy of the aHf[ }\- n O n supercell; E aH f is the energy/atom of 
aH /; E aH f is the energy/atom of aHfO. 
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FIG. 1: Comparison of VASP (large solid circles) and CE (larger open squares, red online) formation 
energies, AEf, and a ground-state analysis on structures with 18 or fewer octahedral-interstitial sites 
(smaller open squares, blue online). Extension of the convex hull towards the formation energy 
of monoclinic Hafnia, HfOz, indicates that the four ordered GS at X=l/6, 1/3, 5/12 and 1/2 
are predicted to be GS of the Hf — O binary. This is not in agreement with experiment which 
suggests that the solubility limit is approximately -ff/Oo.28 (X max ss 0.28; vertical dashed line, 
green online); i.e. in the solid solution Hfi-yOy, Y cs 0.22. 



B. The Cluster Expansion Hamiltonian 

The cluster expansion, CEP% is a compact representation of the configurational total 
energy. In the aHf[ \\-xOx system, the solid solution configuration is described by pseu- 
dospin occupation variables <T;, which take values ai = — 1 when site-i is occupied by [ 
and ctj = +1 when site- 2 is occupied by O. 

The CE parameterizes the configurational energy, per exchangeable cation, as a polyno- 
mial in pseudospin occupation variables: 

E(a) = Y, m ^\H a i) ( 2 ) 

e \ie£' I 

Cluster I is defined as a set of lattice sites. The sum is taken over all clusters I that are 
not symmetrically equivalent in the high-T structure space group, and the average is taken 
over all clusters £' that are symmetrically equivalent to L Coefficients Jg are called effective 
cluster interactions, ECI, and the multiplicity of a cluster, m?, is the number of symmetrically 
equivalent clusters, divided by the number of cation sites. The ECI are obtained by fitting 
a set of VASP FP calculated structure energies, {Estr}- The resulting CE can be improved 
as necessary by increasing the number of clusters P. and/or the number of Estr used in the 
fit. 

Fitting was performed with the Alloy Theoretic Automated Toolkit (ATAT/^EzHSD -which 
automates most of the tasks associated with the construction of a CE Hamiltonian. A 
complete description of the algorithms underlying the code can be found in 28 . The zero- 
and point-cluster values were -0.571537 eV and 0.013973 eV, respectively. The six pair and 
one 3-body ECI are plotted in Figs. [2k and [2b (open symbols, red online), as are ECI for 
ZrOx (solid black symbols) and TiOx (open symbols, blue online). As in ZrOx and TiOx, 
nearest neighbor (nn) — pairs are highly energetic, and therefore strongly avoided; hence 
nn-pair ECI are strongly attractive (ECI >0, for 0-[ ] nn pairs); but beyond nn-pairs, the 
pairwise ECI are smaller; however the 3'rd and 4'th nn pair-ECI in HfOx are significantly 
larger than corresponding terms for ZrOx and TiOx- As in ZrOx, the ratio of ECI parallel 
(J||) and perpendicular (Jj_) to c Hex , respectively, is J\\/J± ~ 2.5; for TiOx, J\\/J± ~ 5. 
These results are similar to those presented by Ruban et al.™ although the ECI presented 
here are not identically comparable owing to different treatments of relaxation energies. 
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FIG. 2: Effective Cluster Interactions (ECI) for pair and 3-body interactions. Solid and dotted 
lines are to guide the eye. Results for the systems TiOx (open blue symbols online), ZrOx 
(solid black symbols) and HfOx (open red symbols online), (a) The first two pair-ECI are for 
nearest-neighbor (nn) 0-[ ] pairs that are parallel- (Jh) and perpendicular (</]_), respectively, to 
the hexagonal c-axis; the second pair parallel to cjjex is J(\ (pairwise-ECI are plotted as functions 
of inter-site separation) , note that the results for HfOx are very similar to those for ZrOx except 
that the 3'rd and 4'th (J/j) nn-pairs are significantly larger in HfOx', (b) 3-body interactions are 
plotted as functions of the index nj_j_fc which increases, nonlinearly, as the area of triangle i-j-k 
increases. Large positive pairwise ECI imply strong pairwise 0-[ ] nn-attractions, i.e. strong 
pairwise O — O nn- repulsions. 

II. RESULTS AND DISCUSSION 
A. Ground-States 



The CE was used for a ground-state (GS) analysis that included all configurations of 
| ] and O in systems of 18 or fewer H /-atoms (octahedral interstitial sites); a total of 
2 18 = 262,144 structures (reduced by symmetry). Five GS were identified in the range, 



< X < 1/2, i.e. at X = 0, 1/6, 1/3, 5/12 and 1/2; solid circles (blue online) on the 
convex hull (solid line) in Fig. [I] The extension of the convex hull towards monoclinic 
hafnia (Hf0 2 ) is also plotted in Fig. [I] The CE-results suggest that all four VASP-GS 
in the aHf[ \\-xOx subsystem are also GS of the H f — O binary. The VASP-predicted 
maximum solubility of O in Hf is X max w 0.5, significantly greater than the experimental 
value of X max rj 0.28. 

Larger open squares (red online) in Figure [I] are CE-calculated values for the AEj that 
correspond to the VASP calculations, and the smaller open squares (blue online) are AEf for 
the remaining 262,144-106=262,038 structures in the GS analysis. All space group determi- 
nations were performed with the FINDSYM program. 125 1 31 1 



TABLE I: Crystal structure parameters for predicted ground- 
state phases in the aHf[ ]i~xOx system. Cell constants 
are given in A. 



System 


X 


Space Group 


Calculated cell 


Idealized 




atomic 


IT number 


constants 


Atomic 




fraction O 


Pearson Symbol 


(A) 


Coordinates 


HkO 


1/6 


R3 


a ps v^flo 


0: 0, 0, 






148 


= 5.5391 


Hf: 1/3, 0, 5/12 




1/7 


hP7 


c«3c = 15.183 


Hf: 0,1/3, 5/12 
Hf: 2/3, 2/3, 5/12 
Hf: 2/3, 0, 7/12 
Hf: 0, 2/3, 7/12 
Hf: 1/3, 1/3, 7/12 


HhO 


1/3 


P31c 


a « y/3a 


O: 1/3, 2/3, 1/4 






163 


= 5.5391 


O: 2/3, 1/3, 3/4 




1/4 


hP16 


cps 2c = 10.122 


O: 0, 0, 

O: 0, 0, 1/2 
Hf: 2/3, 2/3, 7/8 
Hf: 1/3, 0, 7/8 
Hf: 0, 1/3, 7/8 
Hf: 0, 2/3, 5/8 
Hf: 1/3, 1/3, 5/8 
Hf: 2/3, 0, 5/8 
Hf: 0, 1/3, 3/8 
Hf: 2/3 2/3, 3/8 
Hf: 1/3 0, 3/8 
Hf: 1/3 1/3, 1/8 
Hf: 2/3 0, 1/8 
Hf: 0, 2/3, 1/8 



Hf 12 5 


5/12 
5/17 


R3 

148 

hP17 


a « y/3a 

= 10.615 

c p» 6co = 30.366 


0: 0, 0, 1/12 
0: 0, 0, 11/12 
0: 0, 0, 1/3 
0: 0, 0, 2/3 
0: 0, 0, 1/2 
Hf: 2/3, 2/3, 13/24 
Hf: 1/3, 0, 13/24 
Hf: 0, 1/3, 13/24 
ff/: 1/3, 1/3, 11/24 
Hf: 2/3, 0, 11/24 
Hf: 0, 2/3, 11/24 
Hf: 2/3, 2/3, 3/8 
Hf: 1/3, 0, 3/8 
#/: 0, 1/3, 3/8 
Hf: 1/3, 1/3, 5/8 
Hf: 2/3, 0, 5/8 
Hf: 0, 2/3, 5/8 


Hf 2 


1/2 
1/3 


P31m 
162 
hP9 


= 5.5391 
c « c = 5.0610 


0: 0, 0, 

0: 1/3, 2/3, 1/2 
0: 2/3, 1/3, 1/2 
Hf: 1/3, 0, 1/4 
Hf: 0, 1/3, 1/4 
Hf: 2/3, 2/3, 1/4 
ff/: 2/3, 0, 3/4 
Hf: 0, 2/3, 3/4 
Hf: 1/3, 1/3, 3/4 
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FIG. 3: Idealized crystal structures of the four cluster-expansion-predicted suboxide ground-states: 
(a) H/qO; (b) Hf^O; (c) HfaOs; (d) HfoO. Spheres connected by bond-sticks (yellowish-green 
online) represent Hf. Isolated spheres with bond-sticks (blue online) represent oxygen. Isolated 
spheres (red online) represent vacant octahedral sites. 



Ground State crystal structures of the VASP- and CE-GS in Hf — HfO are described in 
Table I and their idealized structures are drawn in Figures [3]a-d: where Hf is represented by 
spheres connected with bond-sticks (yellowish-green online); O is represented by isolated 
spheres with bond-sticks (blue online); and [ ] are represented by isolated spheres (red 
online). As in the ZrOx system, all GS structures are characterized by O — O nn-avoidance 
both parallel- and perpendicular to cuex- 

The VASP-CE-predicted R3 Hf e O GS is the same as the experimental low-T struc- 
ture reported by Hirabayashi et al.™ Space group relations, require a first-order P63HIHIC 
^ R3 disorder '= r order transition between the P63mmc disordered phase and the 
R3 HJqO ordered phase. The Hf§0 GS is the only GS within the experimental solubility 
range < X £, 0.28; all the other computationally predicted GS-phases are presumably 
metastable. 

B. Finite Temperature Calculations 

1. The Phase Diagram 

A first principles phase diagram (FPPD) calculation was performed with conoical- and 
grand canonical Monte Carlo (MC) simulations using the emc2 code which is part of the 
ATAT package^"^. Input parameters for emc2 were: a simulation box with at least 4,050 
octahedral sites; 2000 Monte Carlo passes. The predicted phase diagram is shown in Figure 
111 Most phase boundaries were determined by following order-parameters (77) of the various 
ordered phases as functions of X and T. Dotted boundaries are used to acknowledge uncer- 
tainties in phase boundary determinations. In particular, boundaries of the possible Devil's 
Staircase (DS? in Fig. Ill) regions are are labeled DS?) are poorly defined, and the interior 
structures of these regions are undetermined. 

2. HfoO 

Interstitial ordering of O and [ ] in hep HfOx was studied by Hirabayashi et al. 18 who 
used electron- and neutron diffraction to analyse single crystals with bulk compositions in 
the range of < at% O < 20 (0 < X < 0.25); described as Hf0 1/6 - and Hf0 1/6+ for 
samples with less than or more than one O-atom per six H /-atoms. The structure that 
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FIG. 4: Calculated phase diagram for the system aHf[ ]i_xOx: Approximate regions in which 
the calculations predict a cascade of ordered phases is labeled DS? to indicate possible Devil's 
Staircases of closely related ordered phases. 

Hirabayashi et. al™ report for HfOi/^ has R3 space group symmetry and is identical to 
the VASP-GS at Hf 6 (Fig. [5] and Table I). 

The FPPD-predicted order- disorder transition in Hf§0 (~ 325K) is first-order (Fig. ? ), 
but significantly lower than than the experimental value (~ 700i^ 18 ) or the calculated value 
from Ruban et al. (600Kp2 predicted transition-order not reported). Typically, FPPD cal- 
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FIG. 5: Calculated order-parameter vs temperature (T) curve for the calculated l'st order transi- 
tion in H/qO. 

culations overestimate ] order-disorder transition temperatures so this result is surprising. 



3. Possible Devil's Staircase in HfyO 

The most O-rich structure determination in Hirabayashi et al™ was for a sample with 
X=0.203. The reported structure has P31c symmetry, and is equivalent to the predicted 
HfoO structure [Fig. [3Fb), Table I] except that in the experimental sample, maximal O-site 
occupancy would be ;$ ([ ] .i6, Oo.s-i)- 

Figures |6^i and b are plots of order-parameter vs. T for for Hf^O. The results plotted in 
Fig. [6k, were calculated with the MC-box-size held constant at 4,050 0:[ ]-sites while the 
number of MC-passes is varried. Almost all the order-parameter plateaus are the same for 
different numbers of MC-passes, which reflects the influence of MC-box size on the ordered- 
phase periodicities that are allowed. Note that the transition temperatures from one plateau 
to another are clearly not converged. The results plotted in Fig. |6)d were calculated with a 
constant numbers of MC-passes (12,000) and various MC-box-sizes; i.e. 0:[ ]-sites. Varying 
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FIG. 6: Calculated order-parameter vs temperature (T) curves that evince possible Devil's Stair- 
cases (DS?) of closely related ordered phases at X=l/3 (a,b) and X=l/2 (c,d): (a,c) Monte-Carlo 
(MC) simulations at constant MC-box-size and different numbers of MC-passes at each T, almost 
always find the same set of ordered-phase plateaus (the notation GS4 2000 means the simulation 
was started in the HfyO-GS and run for 2000 MC-passes per T); (b,d) Calculations at a constant 
number of MC-passes, with various MC-box-sizes yield different plateau sequences because box- 
size determines allowed periodicities for ordered phases (the notation GS4 12,000 6,480 means the 
simulation was started in the HfoO-GS and run for 12,000 MC-passes per T, on an MC-box that 
contains 6,480 sites for 0:[ ] mixing). 
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FIG. 7: MC snapshots of ordered phases in the possible Devil's Staircase in HfyO. For clarity, 
only O- (red online) and vacant-sites ([ ]-sites, blue online) are shown: (a) The GS-structure 
with minor disorder. In (b)-(g) there appears to be a competition between two ordering modes: 
1) a striped mode (c,d,f,g) in which layers _L to cnex have single O-rows that alternate with i 
double- [ ]-rows; 2) a triangular mode (b,e) in which layers _l_ to CHex exhibit an ordered array of 
03-nn-equalateral triangles and [ jg-nn-equalataral triangles. 



MC-box-size allows different ordered-phase periodicities; i.e. allows access to more stairs in 
the cascade of ordered phases. 

In Fig. [7^i-g Monte Carlo snapshots are shown at seven different temperatures. For 
clarity, only O- (red online) and [ ]-sites (blue online) are shown. These MC-snapshots 
appear to indicate a competition between two ordering modes: 1) a striped mode (c,d,f,g) 
in which layers _L to cuex exhibit single O-rows that alternate with double- [ ]-rows; 2) a 
triangular mode (b,e) in which layers _L to c# ei . exhibit ordered arrays of 03-nn-equalateral 
triangles and [ ] 6 -nn-equalateral triangles. It is not clear if this cascade of ordered structures 
constitutes a Devil's Staircase, but the results presented in Fig. |6^ and b suggest that it 
does. 

4- Possible Devil's Staircase in HfiO 

In Figures [6J3 and d, one sees similar order-parameter vs. T systematics for Hj^O as one 
finds in Figures^ and b for HfoO; except that the density of plateaus is much greater in 
Hf 2 0, because a 1:1 0:[ ]-ratio allows a greater number different periodic ordered structures 
than a 2:1 0:[ ]-ratio. 

III. CONCLUSIONS 

Ground-State ordered phases are predicted at X=0, 1/6, 1/3, 5/12 and 1/2, but only 
those at X=0 or X=l/6 are likely to be physically realized because the experimental value 
for the maximum solubility of O in hep HfOx is X max rj 0.28. 

Observed ordered phases at X=l/6 and X=0.202P' agree with predicted GS at X=l/6 
and X=l/3 (but with diluted O-site occupancies). 

In the metastable portion of the H f'Ox phase diagram, (0.28 ^ X) cascades of ordered 
phases, possible Devil's Staircases, are predicted for bulk compositions near Hf 3 and 
Hf 2 0. 
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